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ABSTRACT 

A compact formula for the stress tensor inside a self-gravitating, triaxial ellipsoid in an arbi- 
trary rotation state is given. It contains no singularity in the incompressible medium limit. The 
stress tensor and the quality factor model are used to derive a solution for the energy dissipa- 
tion resulting in the damping (short axis mode) or excitation (long axis) of wobbling. In the 
limit of an ellipsoid of revolution, we compare our solution with earlier ones and show that, 
with appropriate corrections, the differences in damping times estimates are much smaller 
than it has been claimed. 
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1 INTRODUCTION 

Most asteroids rotate in the principal, shortest axis mode: their spin 
axes practically coincide with the directions of the maximum mo- 
ment of inertia. Onl y 45 out of almost 5500 entries of the LCDB 
light curve database dWamer et al.l2009l March 2012 version) refer 
to objects that are possible non-principal axis (NPA) rotators, also 
known as 'tumblers' or wobbling objects. With one exception of 
253 Mathilde, tumblers are rather small, with estimated diameters 
below 20 km, but even in this size range they belong to a minority 
among about 2000 objects of this size with known rotation periods. 

Attitude dynamics of asteroids is shaped mainly by gravita- 
tional torques (exerted either systematically by the Sun and gi- 
ant planets, or sporadically during encounters with other bod- 
ies), collisions, optical and thermal radiation recoil torques, i.e. 
the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect, and 
- last but not least - by energy dissipation due to inelastic defor- 
mations. As far as NPA rotation is concerned, collisions and close 
appro aches trigger tumbling (Scheeres et al. 2000; Paolic chi et al.l 
120021') . Small fragments created from collisions of larger objects are 
also expected to start their liv es in a NPA rotation state. The YORP 
effect also excites wobb ling ( lRubincamir2000l :l Vokrouhlickv et al.l 
l2007l:lBreiter et al.l201 ll), whereas - save for possible resonances - 
distant bodies gravitation torques are neutral in this respect. Thus, 
even accoun t ing fo r observational selection effects mentioned by 
IPravec et al.l ( |2005|) . the dissipative damping seems to override 
other effects in most of cases. 

The mechanism of wobble damping was first identified by 
IPrendergast (J958). In NPA rotation, the centrifugal acceleration 
oscillates periodically, deforming each body fragment. The de- 
formation is not perfectly elastic, so some fraction of fluctuating 
strain-stress energy is dissipated during each precession period and 



converted into heat. Draining the elastic energy affects the kinetic 
energy of rotation which also decreases. Thus the rotation axis is 
driven towards the minimum energy state - rotation around the 
principal axis of maximum inertia. The angular momentum, how- 
ever, is not affected by the energy dissipation, as far as we ignore 
thermal radiation and consider the body as an isolated system. Pren- 
dergast provided a general form of energy dissipation rate equa- 
tion for an oblate spheroicQ based upon the solution of 3D elas- 
ticity equations and the assumption that a constant fraction of the 
oscillating part of elastic energy is dissipated at each precession 
period. The latter assumption defines the now commonly adopted 

'Q- model'. 

iBums & Safronovl Jl973h built upon the general idea of Pren- 
dergast using combination of a spheroidal shape for rotation and 
a bent slender beam approximation for elastic energy. Their sim- 
ple estimate of spin axis alignment time is still in use - some- 
times in the version provided by Harris ( 1994). However, some 
scepticism towards it has been brought by observations of as- 
teroids that do rotate around the principal axis in spite of hav- 
ing Burns-Safronov damping time estimate longer than the age 
of the Solar System. Meanwhile, the problem migrated to geo- 
physics (e.g. Chandler wobble damping), rotation dynamics of 
comets and interste ll ar dus t grains physics. The last branch, stem- 
ming from IPurcelll l ll979b . was finally brought back to the dy- 
namics of comets and as teroids with the sequence o f papers by 
Efroimsky and Lazarian I Lazarian & Efroimsky 1999: Efr oimskvl 
,2000: Efroimsky & Lazarian 2000 : Efroimsky 2001 , 2002) . Their 
main point of novelty is the attempt to discuss a triaxial object, 
represented by a rectangular prism (brick), by solving the com- 
ple te, quasi-static stres s tensor equation (Efroimskv 2000). Later 
on, iMolina et al.l ( l2003h issued t he damping model f or a spheroid 
using the same starting point as |Prendergas3 ( Il958[) . i.e. solving 
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equations for displacements. The work of lSliarma et alj ( |2005|) not 
only provides the solution for a spheroid with two different ways of 
estimating the peak elastic energy required for the Q-model, but it 
also offers a long discussion of shortcomings and problems related 
with earlier papers mentioned in this paragraph. 

Trying to combine the YORP effect with a damping mech- 
anism, we first int ended to use the spheroid based model of 
ISharma et al.l j2005|) for arbitrary shap e asteroids. This approach, 
mentioned in lVokrouhlickv et al.l 1200% . became less appealing af- 
ter a closer inspection, because of substantial difference in the dy- 
namics of bodies w ith and without axial symmetry. On the other 
hand, the solution of lEfroimskvl ( I200Q) . albeit referring to a triaxial 
shape, exhibits a number of drawbacks: 

(i) As a consequence of using a non-smooth, brick-shaped ob- 
ject, the solution of stress equations is inexact, with unknown error 
bounds. 

(ii) Compatibility conditions are not fulfilled, i.e. there is no 
displaceme nts field that might produce the strain tensor found by 
Efroimsky jSharma et al]|2005l) . 

(iii) Rotation dynamics is treated by approximate formulae valid 
only in the neighborhood of principal axis. 

Later on, lEfroimskvl J200lh suggested Fourier series involving the 
Jacobi nome as a remedy for the last item, but none of subsequent 
works has implemented this guideline. In these circumstances, we 
have decided to resume the problem at the point where Efroimsky 
has abandoned it, not only using Fourier series to resolve the last 
problem, but also applying the triaxial ellipsoid shape which re- 
solves the first two objections as well. From this poi nt of view, the 
presen t work combines a stress solution in the style of lSharma et al 



( 2005 ) with energy dissipation treatment in the spirit of lEfroimskv 
( 20011) . 



In Section [2] we first formulate the problem of determining 
the stress tensor and enumerate the assumptions, hoping to help 
a reader less familiar with elastic ity problems. Basic facts a re re- 
called accor ding the tex t books of iLandau & LifshitS ( Il959l) . ISaaj 
(l2005l). and 'Wilmanskil ilOld) . Two independent methods (dis- 
placements approach and stress approach) are used to derive and 
cross-check the final expressions of the stress tensor 

The (2-model of energy dissipation is introduced in Section[3] 
and applied in Section |4] to derive an energy dissipation rate for- 
mula. Section [5] presents wobble damping time equations based 
upon the results of Section |4] and some exemplary results. In Sec- 
tion [6] we present the reduction to a specific case of a spheroid 
where a comparison of our solution with those reported in earlier 
works is possible. We use this opportunity to resolve controversies 
concerning drastically different energy dissipation rates in various 
models. 



2 LINEAR ELASTIC MODEL OF ROTATING 
DEFORMABLE ELLIPSOID 

2.1 Basic terms: strain, stress, and body forces 

In an arbitrary reference frame, we consider a body as a dense union 
of material points. Let the set of position vectors r define a refer- 
ence state (configuration); if the points, for any reason, move with 
respect to the reference state, their position vectors will be incre- 
mented by displacement vectors u(r, t), dependent on time t as well 
as on position r, creating a new state with 



The notion of displacements is too general, because it may include 
rigid body motion - translation and rotation of the entire body. The 
rigid body motions are discarded by the introduction of the strain e 
- the dimensionless, tensor quantity describing deformations of an 
infinitesimal volume element in terms of displacements gradient. 

Assumption 1 : The gradient of displacements is small and its 
square can be neglected. 

Under the above assumption, strain tensor e is symmetric by 
its definition 



-[V« + (V«)T], 



(2) 



where the derivatives are taken with respect to the components of 
r. 

Two kinds of forces have to be considered in a continuous 
medium: volumetric forces and surface forces, known also as body 
forces and tractions, respectively. Body forces represent 'external' 
force field. In our case they include self-gravitation and forces of 
inertia. They are defined as a vector field and specified in terms of 
their volumetric density b(r, t) - force divided by the mass of the 
volume element, so the total volumetric force F acting on a body 
with density p{r) is the result of volume integral 

F= \ pbdV. (3) 

Jv 

Surface traction t" is a vector of the force acting on an infinitesimal 
oriented surface, divided by the area. It describes interactions be- 
tween adjacent volume elements or forces applied directly on the 
boundary. The surface is defined by its unit normal vector n, and 
may lie either on a boundary or inside the body. In order to de- 
scribe traction at each possible direction of a plane passing through 
the point r, each component vector of t" in a given basis is pro- 
jected on each component vector of n, creating the Cauchy stress 
tensor T. Thus, traction on a surface defined by n can be obtained 
from T as (unless explicitly stated, repeated index summation as- 
sumed in all formulae) 



'Tjiij. 



(4) 



The units of stress tensor components are those of force per area. 
A good illustration of the two forces nature is a glass of water: 
body force density is constant throughout the volume (homogenous 
gravitational field), whereas tractions define a hydrostatic pressure 
- vanishing on the top surface, reaching maximum at the bottom 
and depending, as a vector, on the direction of the surface element. 

Assumption 2: The deformable body forms an isolated system 
without internal heat sources. 

The consequences of this assumption are numerous. First of 
all, we can use the linear momentum conservation principle in the 
form of Cauchy equation, linking stress tensor, body forces and the 
acceleration of mass particles in an inertial reference frame 



W -J + pb = p{f + u). 



(5) 



r'{t) = r + u{r,t). 



(1) 



If the reference configuration is not fixed in space and still we want 
it to define the reference frame for Cauchy equation, then r should 
be transferred to the body forces b (subtracted) as the density of 
forces of inertia. In case of rotation, some Corolis type terms in- 
volving H may also appear in the right-hand side ( Tokis 1974). 

Assumption 3: Quasistatic approximation. 

Quasistatic approximation results from setting « = (and any 
Coriolis type «) in the right hand sides of Q. T his simplification 
was generally adopted since |Prendergas3 ( Il958h and we proceed 
similarly in the present work. Having included forces of inertia in 
b, so that r = 0, we solve a static equilibrium equation for T 
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V-T = -p6(r,f), 



(6) 



although body forces can be time dependent. The vaUdity of equa- 
tion (|6} can be justified if the solution of the original Cauchy equa- 
tion is a sum of 'free' acoustic waves of high frequency and forced 
vibrations whose frequencies - presumably much lower - come 
from b. Free vibrations can be neglected from two points of view: 
either we consider them to have zero amplitudes at some initial mo- 
ment and then a slow, adiabatic forcing will not excite them consid- 
erably, or - looking forward to the introduction of some dissipation 
mechanism - the high frequency terms will be quickly damped. In 
both cases the stationary regime oscillations derived from l[6]l will 
have amplitudes that differ from the actual ones by a small quantity 
comparable with the ratio of free to forced vibrations periods. 
Assumption 4: Traction-free surface. 

Boundary conditions for the stress tensor will be specified as 
homogenous Neumann conditions on the surface of a deformable 
body 



T/i = 0. 



(7) 



The uniqueness of T as a solution of this boundary problem is guar- 
anteed by general theorems (Saad 2005), but it is not seen imme- 
diately from three scalar equations of Q, even if we add an addi- 
tional property: according to Assumption 2, the angular momentum 
is conserved, so the stress tenso r should be symmetric, i.e. T = T^. 
At this point. Molin a et alj j200 3) felt free to postulate T = 0, which 
was harshly criticized bv lSharm a et al. (2005). 

Assumption 5: Hookean constitutive relations. 

Let us assume that an asteroid is made of a linear, isotropic 
elastic material with adiabatic Lame shear modulus ji and the Pois- 
son ratio v describing the compressibility (incompressible materials 
have the maximum possible v = 0.5). With these assumptions, the 
Hooke's law - serving as a constitutive equation - states a linear 
relation between the strain and stress 

TiJ = 2t^(j^Sijtie + ei,), (8) 
or, conversely, 

., = ^(r,-^.,trT), (9) 
where Sij is the Kronecker delta. 



2.2 Stress approach vs. displacements approach 

Our first goal is to find the symmetric stress tensor T as a solution of 
equations ^ with boundary conditions (O. In the stress approach, 
the problem is solved directly, by assuming some ansatz on T as a 
function of r. But if we accept Assumption 5, additional conditions 
have to be imposed. Strain is a mathematically meaningful quantity 
if there exists a displacement field that generates it through equa- 
tion l[2j. Even without explicit knowledge of u, this is guaranteed 
by Sa int Venant's compatibility conditions ( ISaaj2005l : fwilmanskil 
|20i3) 



V X (V X e) = 0, 



(10) 



providing 6 independent relations between e/j. Through the consti- 
tutive relations ([8}, compatibility equations provide the identities 
that a meaningful stress tensor has to obey in addition to bound- 
ary conditions. In next section we show that (O and l llOt together 
admit a un ique solution for T. The stress approach was applied by 
lEfroimskv. (2000) to the problem of a rotating rectangular prism. 
Yet, the postulated form of T satisfied only the Cauchy equation 



neither boundary conditions, nor compatibility equations could 
be satisfied exactly (the latter wer e not tested at all) an d the level of 
resulting error remains unknown dSharma et al ]l200l . 

Another way of solving Cauchy equation is the displacements 
approach. Using equations ([2} and we convert the first order 
differential equation for stress tensor (|6} into a second degree equa- 
tion for displacements vector field «, obtaining the quasistatic Lame 
or Cauchy-Navier equation 



l-2v 



-pb. 



(11) 



with boundary conditions 

2v(V-«)m + (1-2v) [v« + (Va)'^] n = 0, (12) 

derived from Q. Equations i ll It with only three Neumann bound- 
ary conditions I ll2t admit a solution u(r,t) which is not unique and 
an arbitrary rigid motion may be added to displacements. But since 
the definition of e involves difi'erentiation, the resulting strain and 
stress tensors are uni quely defined regardless of remaining arbi- 
trary terms. Following IPenisov & Novikovl ( Il987i) we will impose 
two special conditions: the volume integral of the displacements 
field should vanish 



X 



udV = 0, 



and the moment of displacements should also vanish, i.e. 



X 



rx«dV = 0. 



(13) 



(14) 



These six conditions aim at suppressing rigid translation and ro- 
tation terms in displacements and allow a unique determination 
of «, which is of minor interest for the stress tensor recovered 
through l|2j and l[8j, but gives more insight into the question of 
the reference configuration choice and simplify energy and mo- 
mentum balance discussion. Up to the ambiguity in the last tvyp 
conditions, the displacements approach w as taken by |Chree(1895^ 

ilsi 



iDenisov & Novikovl ( ll987h . lMolinaetalJ ( l2003h . and lsharma et al] 
( I2OO5I) . 



We can also observe that restoring the term pii in equa- 
tion i ll It . we obtain a quantitative measure of quasistatic approx- 
imation error. The homogenous solution will involve a frequency 
close to 



Wf 



(15) 



where a is the radius of an object. If body forces are periodic with 
frequency CI (the precession frequency in our case) then, with /i of 
the order of 10 GPa, the ratio Cl/cof may be safely considered small. 

2.3 EUipsoid stress solution 

2.3.1 Homogenous ellipsoid body forces 

Let the reference configuration be a homogeneous rigid ellipsoid 
with semiaxes c Kb ^a. Its shape will be described by two dimen- 
sionless parameters 

b c 

hi = -, 112 = -, (16) 
a 

both taking values < /?,■ < 1. In the reference frame whose cen- 
tre coincides with the centre of mass and the basis vectors e; are 
directed along the principal axes, the parametric equation of the 
interior reads 
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r = qa (sini?cos(/iei sinj?sin0e2 +hi h2 cos&ey,), (17) 

where 0<g<l,0<)?<:nr, 0<^< 2n. The boundary is specified 
by ^ = 1 and the unit normal vector on the boundary is 

n = 0(/! i,/z2, (sin i?cos0ei 

sin)?sin(^e2 +'^7'''2' cosi?e3), (18) 

where <I) is some nonzero function; its explicit definition is not re- 
quired for the vanishing traction condition JT). According to the 
postulate (13}, the centre of ellipsoid remains the centre of mass 
even in a deformed state. 

Given an arbitrary function F(r), the volume integration rule 
for an ellipsoid is 

FdV = a^h\h2 Aq siniJdi? Fd^, (19) 

V Jo Jo Jo 

where F(r) should be expressed in terms of q,&,^ according to 
equation l ll7t . 

The body forces acting on a freely rotating ellipsoid include 
the forces of inertia due to rotation, with force density vector 



bin = rx<b + cjx (rxo), 



(20) 



where the time derivative of rotation vector co is given by Euler 
equations 



l-hi 



l-hi 



(22) 



LL,-~- — -^a)2mei + - — -j-^aji(D3e2- - — -7^11^263, (21) 
1+^2 l+Aj/jj l + /lj 

and the moments of inertia /, for an ellipsoid with mass m 

/i ^1(^+^2) h i + hjhj /3 l + /if 
ma^ 5 ma^ 5 ina^ 5 

are substituted. 

Gravitation inside the ellipsoid results in body forces density 
*gr = -n^fl -r2>'e2- 732^3, (23) 
with constants 

72 = ^RAhhlhjhlh^), (24) 

73 = ^RAhh^i,hlhl,hlhl), 

expressed in terms of gravitation constant 7 and Carlson's elliptic 
integral 



R j{u,v,w,p) : 



ds 



2 Jo (p + s) V(m + s)(v + s)(w + s) 



(25) 



The body forces are linear in coordinates of reference config- 
uration, so we write them as 



b = bin + bgi = Br, 



(26) 



with coordinates independent matrix B having elements 



2 2 
(^2-1-013-71, 



S22 = <^l + '^\-72, 
S33 = W[-fW2-73, 

2a)i(Li2 

B\2 = -J, 

1 

2(Di(jJ3 



(27) 



B23 = - 



2(jJ2<i'3 



and 



-621 



S31 



B32 



■■hlB2i. 



(28) 



2.3.2 More assumptions 



Assumption 6: Displacements and their partial derivatives are small 
quantities of the first order. 

This stronger variant of Assumption 1 allows to treat the prob- 
lem in terms of a first order approximation. Namely, considering 
boundary conditions we can impose them on the reference ellipsoid 
surface, neglecting the displacements that deform it. It also means 
we do not restrict reference ellipsoids t o a figure of equilibr ium 
type solutions like e.g. Jacobi ellipsoids dChandrasekhadl 19691) . In 
the same spirit, we ignore the variations of density due to strain and 
use a constant, mean p whenever it serves to define displacements 
or stress ~ either explicitly, or indirectly (like in body forces ftg,). 

It turns out that postulates l |13t and I ll4t are inherently re- 
lated with the choice of reference configuration that satisfies our 
assumption. In the first approximation, i.e. evaluating volume inte- 
grals over the homogeneous reference ellipsoid, we interpret l |13l l 
as a postulate that the centre of mass position is not altered by dis- 
placements. Similarly, equation I ll4t indirectly leads to the state- 
ment, that displacements velocities do not contribute to the angu- 
lar momentum of the system. Observing that « will depend on Bij 
in exactly the same form as u depends on Bjj (the Cauchy-Navier 
equation is linear and does not involve time derivative), we find 



X 



p(rxu)dV = 0, 



(29) 



as a consequence of l ll4t . provided p is constant and the same 
bounding surface is used in both integrals, i.e. within the first or- 
der approximation. In other wo rds, the postulate l|14t imph es that 
we use Tisserand's mean axes jMunk & MacDonaIdlll96(ll) as the 
reference frame and they approximately coincide with the princi- 
pal axes of the reference ellipsoid. Such choice has a property of 
mini mizing the displacem ents and their velocities. 

ISharma et all ( |2005|) postulated a pre-stressed state of their 
reference spheroid and dropped the constant part of the stress due to 
self-gravitation on the onset of their derivation. They were not con- 
sequent in this point, because they did not do the same with a mean 
part of centrifugal stress. For typical asteroids both effects may be 
comparable. They may even mutually cancel. Thus we do not find 
the pre-stressed state assumption necessary for asteroids, although 
it is important for major objects, like the Earth, where it was origi- 
nally introduced (Love 1934). On the other hand, its role would be 
to provide the rationale for the validity of the six assumptions we 
have already made. 



2.3.3 Displacements approach 

The problem of finding displacements and stress tensor for a 
freely rotat ing ellipsoid (without self-gravitation) was solved by 
lOenisov Novikov (1987). Their impressive solution, apparently 
worked out without the support of computer algebra, was based 
upon the displacements approach, which helped in establishing the 
influence of disp l acem ents field on the moments of inertia. Earlier 
results of IChred ( Il895 ) included the self-gravitation, but only the 
principal axis rotation was considered. We used the stress tensor of 
lOenisov & Novikovl ( Il987h as a test of our solution. 
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The main advantage of using an ellipsoid is the possibility of 
finding the displacements field that satisfies equations i ll It . l ll2t 
and l |26| l as a sum of two homogeneous polynomials of degrees 1 
and 3 i.e. a sum of 13 monomials with dimensionless vector coeffi- 
cients / 

+/210x2v + /30»x3), (30) 

which involves 39 arbitrary constants: 9 for the linear part and 30 
for the cubic terms. It is instructive to add B2\, S31, and S32 as ad- 
ditional, unspecified parameters, raising the number of unknowns 
to 42. Our choice for the form of u is straightforward, although 
not necessarily optimal; for example. lSharma et al. I diool used the 
spherical harmonics basis and reported a smaller number of unde- 
termined coefficients for a spheroid and in an unpublished solution 
for the ellipsoid (Sharma, private communication). 

Using the 'brute force' attack, allowed by the use of an alge- 
braic manipulator, we formulate 42 independent conditions involv- 
ing linearly the coefficients of u and unspecified elements Bij of 
matrix B - the ones with i < j are parameters, and those with i > j 
are unknowns. First, we observe that the centre of mass condition 
il3\ is satisfied identically by i30l . Actually, it has been used im- 
plicitly to drop coordinate independent and quadratic parts of u. 

(i) Equation dl It generates three equations linear in x, y, and z- 
This amounts to 9 equations linking each Bij with 5 coefficients of 
the cubic part of u. 

(ii) Postulating equation we obtain 3 equations - each link- 
ing two coefficients of the linear and six of the cubic part of u. 

(iii) The remaining 30 equations result from boundary con- 
ditions 1 II2I 1. They can be derived either directly, i.e. from 
the polynomial form of ellipsoid surface and normal vector 
jPenisov & Novikovll 1987b . or by considering trigonometric poly- 
nomials of (f) and & resulting from the substitution of equations UH 
and l llSt . In the latter case, we first equate to the coefficients 
of sin3(/i, cos 30, sin 20 and cos 20 which have a common factor 
(sinj?)^ or (cosi?)^; this provides 12 conditions. Then, in the co- 
efficients of sin0 and cos0, we separate the terms with sm§ and 
sin 3)?, and set them to 0, obtaining another 12 conditions. Finally, 
in the part independent on (p we separate terms factored by cost? 
and cos3i? resulting in the last 6 conditions. 

Solving the Cramer system of 42 linear equations we obtain a 
unique solution for the 39 coefficients of u as linear combinations 
of Bi j with I < j. For the remaining three unknowns we recover 
equation l |28t which becomes the condition of existence of the dis- 
placements solution OOt for a n ell ipsoid - slightly more general 
than the condition found bv lPenisov & Novikovi ( Il987l) . 

Unfortunately, the resulting expressions of displacements field 
H are too long to be quoted. Each Bij appearing in / has its own 
multiplier - a rational function of h^, and v. The only common 
factor of all / vectors is pa^ii~^ = wjT^. This indicates a link be- 
tween Assumptions 3 and 6. 

2.3.4 Stress approach 

The structure of displacements field is inherited by the stress 
and strain tensors. Each Tjj and eij is a linear combination of 
Bij, including zero and second degree monomials of x, y, and 
z. We focus the discussion on the stress tensor, because T can 
be of interest in studying the breakup of spinning asteroids (e.g. 



IWashabaugh & Scheeresll2002l) ). whereas strain components eij in 
the elastic material are easily found from T using constitutive rela- 
tions l|9j. 

Introducing 

hi hih2 
to benefit from some symmetries, we can write the general form of 
T for the elastic ellipsoid as 

-xyAi2-yzA23-xzAi3), (32) 

where pa^A represents the stress tensor at the origin x = y = z = 0, 
and all matrices are symmetric. Thus the stress definition requires 
42 matrix elements to be determined. Their expressions resulting 
from direct substitution of displacements solution are unwieldy, so 
we decided to apply the second possible way of determining T - 
the stress approach. Quasistatic equilibrium condition ^ does not 
involve A and leads to 9 linear relations between A'-' and body force 
matrix B. Boundary conditions Q generate 30 relations between 
the elements of A and A'^. However, the resulting set of 39 linear 
equations admits solutions only if conditions l|28) are satisfied and 
then its rank drops to 36 allowing a unique solution for all A'-' as 
linear combinations of A and B elements. 

'Central stress' A results from compatibility equations MQ\ , 
forming two independent subsystems of equations: 





d^ei2 


3^623 


5^613 


dydz 


dxdz 


5x2 


dxdy ' 


d^e22 


d^e2-i 


3^613 

5y2 


d^ei2 


dxdz 


dxdy 


dydz ' 


3^633 


d^en 


a^ei2 


5^ £23 


dxdy 


dydz 


&2 


dxdz ' 


and 








23^<?l2 


d^en 


3^622 
^ 5x2 ' 




dxdy 






2'9^e23 


d^e22 


5^633 
5>'2 ' 




dydz 


dz^ 




2<9^ei3 


d^en 


3^633 
5x2 ■ 




dxdz 


dz^ 





(33) 



(34) 



Relating e with T by means of (|9}, substituting the general form 
of stress l l32t and making use of A'-' expressed in terms of A and 
B, we find that subsequent equations of ( 1331 ) directly define A23, 
Ai3, and A12, respectively, in terms of S,j with the same subscripts. 
On the other hand, equations l l34b form a system of coupled linear 
equations for A,-,- with right hand sides depending on B,-,-. Its solution 
takes form 



( 2Au ] 




' fill - 


2hfA22 


= (i+L-1r) 


B22 


. 2/!j2A33 , 




[ S33 j 



requiring the inverse of a 3 x 3 matrix L. Although L is not compli- 
cated and invertible by elementary means, an explicit formula for 
L^' is too long to be explicitly quoted. The detailed solution for A'-' 
and matrices L, R are given in Appendix IaI 

We have verified that both methods (displacements approach 
and stress equations with compatibility conditions) lead to the 
same final results. Nevertheless, the reduction of T to the simple 
form given in Appendix|A]starting from the displacements solution 
would be a tedious exercise. 
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3 THE Q MODEL OF DISSIPATION 



3.2 Quality factor principle 



3.1 Energy of a deformable ellipsoid 

For a rigid ellipsoid in free rotation with angular velocity co, the 
only part of energy that matters is kinetic energy 



If 9 1 

^0=r piojxr)-dV = -ojilijoij. 



(36) 



conserved during the motion. In the principal axes system, the ten- 
sor of inertia I is diagonal, i.e /,j = Sijli^ given in equation J22b . In 
the absence of external torques, the angular momentum 



H= I prx(ojxr)dV = lijcojei, 

JV 



(37) 



is also conserved in an inertial frame, i.e. provided we account for 
the rotation of body-fixed basis vectors e,. Only the norm H = \\H\\ 
is constant in the body frame. 

In a deformable body, the velocity of displacements « adds 
up to the total kinetic energy K and angular momentum. Moreover, 
the results of volume integration are affected by the fluctuations of 
density and by the bounding surface deformations. However, under 
Assumption 2, as long as the deformations are elastic and the cen- 
tre of mass is not altered by displacements, the angular momentum 
remains constant in the fixed frame. Total kinetic energy may vary, 
but the sum K + U, where U is the potential energy of deforma- 
tion, remains constant. A rigorous discussion of energy exchange 
between K and U can be found in lMunk & MacDonalj i ll 9601) or 
iLam beck ( 1980). For the present discussion, we approximate the 
kinetic energy as K x Kq and consider potential energy U ~ Us, 
storing the work of both body forces and tractions, to be the elastic 
energy 



edV^O, 



with the energy density e defined as 



1 



(38) 



(39) 



which is evaluated, according to the first order approximation, with 
the volume integral taken over the reference ellipsoid. Under the 
Assumption 5, we can also expre ss e in terms of the stress tensor 
alone, obtaining l lEfroimskvll2000l) 



1 



y(trT)2 
1 + v 



(40) 



Inelasticity disturbs the ideal picture of Hookean oscillations. 
Diff'erent rheological models have been constructed in hope to ad- 
just constitutive relations between stress and strain to the reality 
of the material world. Extended constitutive relations are formu- 
lated either as differential equations, involving T and/or e, or in 
terms of integrals of creep and relaxation functions. This leads to 
the occurrence of a time lag between forced oscillations of strain 
and stress. In this respect, the situation becomes similar to a pe- 
riodically driven harmonic oscillator with damping: due to a lag 
between velocity and position, the time derivative of potential en- 
ergy integrated over a forcing cycle does not vanish and generates 
a power deficit credited by the external forcing in order to sustain 
stationary oscillations. In the case of our study, the power supply 
comes from the kinetic energy A^o and is not unlimited. Moreover, 
a coupling exists between the power supply and demand, because 
the amplitude of stationary vibrations depends on the excess of ki- 
netic energy over the ground state of 5/301^. 



Introducing the quality factor Q as an empirical parameter, one 
may, in principle, discuss the energy dissipation in vibrating in- 
elastic materials without any explicit knowledge of their consti- 
tutive relations. In practice, however, the problem of unknown 
rheology leaks into the question of a proper definition of Q 
and of its dependence on driving frequency, temperature, etc. 
( Efroimsky & Williams 2009). In the limit, a perfect definition of Q 
is probably not easier than defining an adequate rheological model 
and solving the related inelast i c osci llations problem. 

lO'Connell & Budianskvl ( Il978h tried to put some order into 
a growing number of different Q definitions. They warned about 
'confusion between some ill-defined 2 of a process' and 'an in- 
trinsic Q of the material'. Their generally acclaimed definition of 
quality factor as the ratio of real and imaginary parts of compliance 
leads 'for a large class of viscoelastic materials' to the common 
rephrasing as 



e = 2;r 



(2£av) 



(41) 



where AE is the energy lost during one period of a harmonic (pure 
sine) loading, and £av is the 'av erage stored energy' during the 
loading cycle jO'Connell & Bud ianskv 1978). The difficult point 
of this apparently simple definition is how to interpret 

tav ttt some 

particular problem, even if we use the strain energy I l38l39t for this 
purp ose. 

Isharma et aD (l2005h complained that 'definitions of measures 
of energy fluctuations corresponding to the type of loadings en- 
countered with tumbling bodies are not readily available', meaning 
the presence of a constant term in body forces that are not purely 
periodic. Although their main solution followed the usual habit of 
dropping the constant part from strain and stress when plugging 
Us into the definition l l41t . they express serious doubts and find 
'no easily identifiable reason' for it, except that of comparison with 
earlier works. They proposed an alternative approach including the 
effect of the avera ge (b). When starting the present work, we did 
share the doubts of Sharma et ai] ( I2OO5I) and were attracted by the 
alternative way of estimating the fluctuating energy amount. But, 
having rejected them in later stage, we feel obliged to explain our 
point of view. 

Consider the example given bv lSharma et al.l ( l2005h : a scalar 
stress T = a^ + ai sinf. There is no doubt that (t^^ with # is dif- 
ferent than for ao = 0- Moreover, with tO the plot of L'e(f) can 
be dominated by the sin t, whereas dropping the constant part we 
obtain only cos 2t. However, a similar situation is met in a damped 
harmonic oscillator driven by + ^1 cosf: the potential energy of 
stationary solution has the same dependence on ao, but the power 
loss is completely independent on aq- The statement that only a 
time variable part of st ress may dissipate t he energy can be found 
already in the paper of IPrendergastI ( Il958l) . There is no reason to 
doubt it. Thus the only question is: does the the dissipation by pe- 
riodic part of the stress depend on the constant part ? We cannot 
rule out such possibility if relations between stress and strain are 
strongly nonlinear or the dissipation mechanism is complicated, but 
- on the other hand - it is quite unlikely that in these circumstances 
a simple term apai s inf will properly describe the dependence, as 
ISharma et iD l l2005l) suggest, mentioning a friction due to grain 
boundary or crack surface sliding. So, the alternative estimate of the 
stored energy, as proposed by Sharma, neither looks promising far 
from almost-linear, weakly damped elasticity model, nor it behaves 
properly within this area - the latter readily seen if we consider the 
harmonic oscillator. Moreover, it does not account for the fact, that 
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a k-\h harmonic in a general harmonic load works k times during 
the fundamental period - the comment that ap plies both the main 
and alternative solution o f ISharma et id] ( l2005h . As a minor remark 
we can add the contradiction between qualifying gravitation as an 
ignorable pre-stress and mean centrifugal force as a factor that con - 
tributes to the energy dissipation, present in Sharm a et al. 
The alternative recipe for E^^ has a mathematical meaning, but in 
our opinion, it will define an alternative '2 of a process' which may 
be too far from the 'intrinsic Q" in terms of numerical values and 
dependence or independence on parameters. 

We believe that the main line of rea soning and the warnings is- 
sued by jO'Connell & Budianskvl[r978h are sufficient to formulate 
the recipe for £av leading to a reasonable, material-based quality 
factor. The rule seems to be: stay close to the property that, for 
sufficiently high Q values. 



■ tani/), 



(42) 



where (p is the phase lag angle between the stress and the strain 
(so-called loss angle). This rule validates the separation of contri- 
butions from subsequent harmonics of body forces b, as well as the 
r ejection of it s const ant part in the stress solution, i.e. the procedure 
of lEfroimskvl (l2000'). Note that, for a single harmonic, taking twice 
the mean value of its square we obtain its squared amplitude which 
explains why 2£av in equation (141b is often replaced by 'peak en- 
ergy' Bp dO'Connell & Budianskvlll978l : lLambecklll98(]h . 

Energy is dissipated by the work of body forces. The basic 
formula for the work rate W reads 



W 



. r edv= r 

JV JV 



TijeijdV. 



(43) 



In the conservative, elastic case, when stress and strain are defined 
by equations l l32t and depending on time-periodic functions 
Bjj, the integral over the fundamental period of wobbling renders 
no work, because each term of the sum Tijejj is a purely periodic 
function of time. 

Assumption 7: Inelastic oscillations will be described by the 
following, heuristic approximation: 

• Periodic terms of Bjj in stress tensor T are taken directly from 
the definition J27t and J28t . 

• Periodic terms of 6,j in strain tensor e, derived using equa- 
tion l|9), are modified by adding a phase lag (f to the argument of 
each harmonic, and dividing the amplitudes by cos^. 

• The phase lag is independent on coordinates x, y, z, and related 
with the quality factor by equation i42l . 

• The quality factor is independent on the frequency of forcing 
terms. 

Let the fundamental frequency of wobbling be CI, with an associ- 
ated period 



27T 

n"' 



(44) 



and consider an exemplary product of periodic terms appearing in 
equation l l43t under Assumption 7 

/ \ d / 

Pex = {cpCos(pQ.t) + SpSin{piit)) — cos {qQt-w) 

^ ' dt\ cos (fi 



H sin (qQ.t -<fi) 

cos (fi 



(45) 



Straightforward computation leads to the conclusion, that the time 
integral over the period P vanishes for p i= q, whereas p = q t 



leads to 



Pex dt = (cp + Jp) tan <^ = 



= p^ (^{c pcos {pnt) + sp sin (pClt)fy (46) 

This example establishes the link between Assumption 7 and the 
operational rule of computing the energy loss due to the work per 
cycle P as the sum 



A£ = - rVdf = -|2p(2{/,) 



(47) 



in agreement with equation i4U . where Up is the part of elastic 
energy Ue from equation ( I38t involving only the p-th harmonic in 
each Bij term entering strain and stress. 

An recent in-depth critica l review of some issues concerning 
Assumption 7 can be found in l lEfroimskvll2012l) . 



4 ENERGY DISSIPATION RATE 
4.1 Volume integration 

Energy dissipation rate can be obtained from equation ( 147 b as 
AE 2Q. 



^ P>i 



(48) 



Having assumed the independence of ip and Q on coordinates, 
we can perform the volume integration required for Up before 
the time average. This allows a considerable economy of expres- 
sions; substituting the general form of T from equation i32l into 
equation ( l4Qt . making use of the expressions for A, A'-' from Ap- 
pendixjA] plugging in the definition of B in ( |27|28t , and integrating 
over the ellipsoid volume according to equation ( |19t , we find 



i^p) = ^(-n(K£)-*-«22(N£)-*-33(H 

+°'23{[<^l\[oji\)+Pl2{loJl<^2fp) 

+13x3 {Wic^if,) +P23 {lc020J3]l)) , (49) 

where m is the ellipsoid mass and aij, Pij are dimensionless ratio- 
nal functions of hi, hi and v. For any function F represented as 
Fourier series, the symbol \E\p designates its p-th harmonic, i.e. a 
trigonometric monomial with argument pQx. 

Regretfully, the full form of ct/j and j8,j is too long to be explic- 
itly quoted (although short enough to be efficiently programmed), 
but the expressions are available from the authors in an electronic 
form. 

4.2 Fourier harmonics of angular velicity 

Equation ( |49t requires the recall of basic facts about the free rota- 
tion of rigid body ( Whittaker 1952) for the specific case of a ho- 
mogenous ellipsoid. As usually, we have to distinguish the Short 
Axis Mode (SAM) of rotation, when angular velocity vector o) cir- 
culates around 63, and the Long Axis Mode (LAM), when <u circu- 
lates around ei . Quantities referring to the SAM will have subscript 
i = 3, and those for the LAM - subscript i- = I (in this paper, we 
use s exclusively for labeling the rotation mode). 
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Using two invariants of the free top, i.e kinetic energy A^o and 
angular momentum H (equations l l36t and ll37t). with diagonal ma- 
trix of inertia I given by equation p2ll. we define an auxiliary quan 
tity 

2Ko 



: mertia I given by equation \22\. we define £ 
[ dDeprit & ElipJll993l:lBreiter et alj|201lh 



(50) 



(51) 



such that (33 < ^ < Qi , where 
ai = i;\ 

A nominal angular velocity of rotation Haj is often adopted as 
a scaling factor in the energy dissipation models. Its value would be 
equal to co only in the principal axis rotation around ej. We prefer 
to use a quantity appropriate for the LAM as well, so let us define 
the nominal angular rates as 

(bj = Hcij, (52) 

such that each ajj = coin the principal axis rotation around ej. 
Let us define 



«3 = ^/(a\ -y{){a2-a-i), m = V("i -"2)(-^-"3)- (53) 

Then, the components of angular velocity cl), in the body fixed 
frame are expressible in terms of the Jacobi elliptic functions with 
argument t, = Hn,, t and modulus fcj, where 

ki = — = — 

Thus, in any rotation mode. 



(54) 



CJi = toll 



0)2 = (JJl 



£013 = tLi3 



r 


-03 




— (33 


t 


-fl3 




— (33 






J ai 


— 03 



' S 



(55) 



with specific functions 

f = ±cn(T3,/t3), Ff' = sn(T3,A:3), F*^' = ±dn(T3,yt3), (56) 
and 

= ±dn(Ti,yti), =+fcisn(ri,yti), Ff *= cn(Ti,yti), (57) 

for the SAM and LAM, respectively. 

The fundamental frequency of wobbling, appearing in equa- 
tion i44[ as n, IS 

nHn^ 

n. = — , (58) 

where stands for the complete elliptic integral of the first kind 
dx 



K, = Kiks) - 



yjl-k^ sin^ X 



(59) 



Similarly, we will use = E{ks) for the complete elliptic integral 
of the second kind. 

In the present work, we are interested only in the squares and 
products of angular velocity components. Given the Jacobi elliptic 
functions with argument r., and modulus ks, we can expand their 
squares and products in Fourier series of angle 



(60) 



From the expressions available in iBvrd & FriedmanI ( I1954I) , we can 
easily derive series 



^(Ts,ks) = XQ-^X2pCOs2ptl/s 



^(T,,A;,) = l-Xo + J^X2pCos2ptff„ 



(61) 
(62) 



dn 



=(T„^,) = I-fc2Xo + fc2^X2pCos2piA„ (63) 



sn(Ts,ks)cn(Ts,ks) = ^Y2pS'm2pi//s, 



p=l 



cn(Ts,ks)dn(Ts,ks) = /rj^X2p_i cos(2p- l)i/'s 

p=l 

CO 

sn(Ts,ks)dn(Ts,ks) = ks^Y2p-is'm(2p- l)if/s, 
p=l 



4 1-^ 



with coefficients 
^0 



k^\ K, 



\2 . j/2 



i + qi 

involving the Jacobi's nome 
q^ = exp(-;rK;5/Kj), 



(64) 
(65) 
(66) 

(67) 
(68) 
(69) 

(70) 



where K'j = K(k'.,) is the elliptic integral with complementary mod- 
ulus 



(71) 



Quickly convergent series for the nome 

q, = ^,(1 +2^ + I5d + 150^^2 + 1707 f]*- + ...), (72) 
with 

1- 



(73) 



2(1+ Vfel)' 

can be used dinned [l902l : iBvrd & Friedmaiilll954h . Even close to 
the separatrix x 02), the first five terms provide a relative error 
of 10"^ for ks = 0.999. However, for the values of ks close to 1, it 
is better to use 

1- 



2{l+ ^/k;y 



(74) 



in the right-hand side of ( 1721 ) instead of and obtaining the com- 
plementary nome q^, which may serve to compute q, through the 
relation 



Inq j Inq, = n^. 



(75) 



An additional benefit of the nome is also a quickly convergent se- 
ries for the elliptic integral 

l + 2j]qf . (76) 

P=l . 
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The mean values required in equation l |49t follow directly 
from the presented Fourier series. For example, in the SAM 



2p 



01-03 



((^2;,COs2pi/'3)" 



(I 



and in the LAM 



2 \ai-aij 2P' 
0, 



2]2 \_^\ / -yi-«3 \' ,,4y2 



(I' 



\2p-\ 



(77) 
(78) 

(79) 



with X2p depending on q3 and qi , respectively. Note, that some of 
the mean values can be negative, like 



(N]2.N]2. 



-2-2 



2 (ai-aT,){a2-a-i) ^P' 



(80) 



in the short axis mode. However, their associated 0-,^ are also nega- 
tive, so there is no subtraction in equation l |49l l. 



4.3 Final expressions 

Performing the necessary substitutions in equation ( I48t . we find an 
expression for the energy loss rate 



4 -5 



with dimensionless 



(81) 



T3 



"Pi 



(82) 



Z\ {Plih )Mi3 + P2(h )Mn + Piih )Mo + PAih )M23) , 



(83) 



(84) 



(note the swapped M23 and M12), where 

Q 5 nil, 
Zs = — = —. 

His 2a.5Kj 

First, we recall that the leading factor depends on the semi-major 
axis of ellipsoid a, its mass m, density p, the fifth power of the nom- 
inal rotation rate Cbs (resulting from the division of angular momen- 
tum H by the related moment of inertia), on Lame shear modulus yu, 
and quality factor Q. Functions Pi(ks) depend on the ratio of kinetic 
energy and angular momentum through an elliptic modulus that 
enters the Jacobi's nome qj, and have the form of infinite sums 



Pl(k.) 



P2(k.) = 



t» ,,3 2n-l 

y (2p-l)-^q/ 



2;a2 ' 



Piiks) 



;^(l-q?'0 



i2p 



(85) 



(86) 



(87) 



(88) 



although in practice only a few leading terms should be sufficient. 
Finally, Mjj and Mq are dimensionless, positive coefficients de- 
pending only on the shape, (through hi, /12) and on the Poisson's 



ratio V. In terms of the coefficients from equation ( |49t , they are 



Mo 



d\2d\3d2idij ' 



16 



^12^^13^23 



a,rf23Q'll 



«Vl3Q'22 a^dna^i 
+ — + ■ 



2 2 2 2 

d\2 di3 



dndn ^12^23 



2 2 
fl2'J3Q'23 



dnd23 



d23 



(89) 



(90) 



where dij = (a, - aj). Appendix IbI contains the full expressions of 
Mij and Mq. For the reasons explained in the next section, we give 
them with the fixed Poisson's ratio v = 0.25. 



4.4 Poisson's ratio 

All recent models of spin axis relaxation assume the Poisson's ra- 
tio V = 0.25, i.e. equal Lame constants A= /i. Authors justify it by 
the fac t, that this i s app roximately a typical value for most of cold 
solids ( lEfroimskvl[200Q) . Earlier, IPrendergas^ lll958[) con s idered 
an incompressible object with v - 0.5. Onlv lMolina et alj j2003h 
maintain the explicit dependence on v in their final formulae for a 
spheroid. 

The present model also maintains v in the final expressions, 
so we are in a favorable situation to estimate the sen sitivity of Es 
on its value. Interestingly, in contrast to the results of lMolina et al.l 
12003), the dependence of Mjj and Mq on the Poisson's occurs to be 
very weak. As a function of < v < 0.5, the values of M coefficients 
vary o n the level of at mo st 10"^ (relatively), whereas the solu- 
tion of lMolina etai] ( |2003|) exhibits the dependence on the level of 
10-1. This property came unexpected, because a, ; still contained a 

9—1 

factor (1 - y ) , that later vanished in Mq. In these circumstances, 
we fix the value of y = 0.25 as a physically realistic one, which 
considerably simplifies expressions, but the results will fairly well 
apply to an incompressible case with y = 0.5. 



5 WOBBLE DAMPING TIME 

Let us define a 'wobbling angle' Gs as the maximum angle between 
the angular momentum vector H and a relevant axis {Oz in SAM, 
or Ox in LAM) attained during the wobbling cycle of a rigid body, 
namely 



max arccos 



H 



(91) 



In the principal axis rotation = 0, regardless of prograde or ret- 
rograde case. The other limit, Os = 90°, refers to unstable rotation 
around an intermediate axis of inertia Oy, or to the nonperiodic ro- 
tation on the separatrix, which falls beyond our model. 
There is a direct relation between and the variable 



Jls = a2- {02 -as) cos, 



(92) 



We have introduced the subscript s Xo Ji m order to facilitate the 
distinction of modes, although the primary definition l |50l l is univer- 
sal. The modulus ks is also related to the wobbling angle through 



^3 



sin 05 


Vi 


+ Ks cos2 9 s 


1 


02 - 




ai-a2 



1-^ 



(93) 
(94) 
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Since the dissipation of energy does not affect angular mo- 
mentum, and the energy is drained from the kinetic Kq, the differ- 
entiation of equation \50i leads to 



2Es 



2a p mibi i 



(95) 



where Es has been taken from equation l lSlt . 
On the other hand, 

= 2 ia2- as)sm6sCosds9s, (96) 

and we can equate the two relations, obtaining the differential equa- 
tion 

de, E, 

df {a2 - a J ) sin 0s cos d,, 

The resulting quadrature gives the time Ts required to change wob- 
bling angle from the initial 9j to the final 6'^ 



(97) 



(as-a2)fiQ r^' sin 6/,; cos 6*. 



(98) 



where "P, should be expressed in terms of the wobbling angle using 
equation i93\ . 
In particular. 



73 = - 



7 -3 



/,2(l + /,2)(l-/,2) 



5{l+h^,hl) 



r 



sin 03 cos 03 
^3 



d03. 



implying 6'^ < S^, and 



7 - 3 



;,2(l-/.2)(l+/.2) 



5{\+h\hfi 



sin^i cosf? 



de. 



(99) 



(100) 



with 



'1 < ^'l 



For the reference with earlier works, we will use a 



shape parameter jSharma et alj2003) defined as 



D,(hi,h2) = T, 



2 - 3 
a pu)^ 



(101) 



for prescribed integration limits. 

Of course, the notion of wobble damping time is properly 
related to or, if the evolution starts in the LAM and continue 
through the SAM, to the sum Ti +72- In the long axis mode, en- 
ergy dissipation excites wobbling, driving the angular momentum 
vector towards the separatrix. A term 'excitation time' seems more 
appropriate for Ti . 

As an illustrative example, we first plot damping/excitation 
times for a family of ellipsoids with hi = h2 = h, assuming sample 
physical data a = 1 km, p = 2000 kgm-\ = 10^ Pa, and 2 = 100, 
chosen more for an ease of scaling than for the reference to some 
specific case. We have assumed cii^ = In/W^ for and computed 
an equivalent 



l + Zi^ 



(102) 



to be used in Ti . 

Interested in a possibly complete history, we start the evolu- 
tion at the LAM, with 6*'' = 5° and integrate equation JlOOt up to 
0j = 85°. The results are displayed in Fig.[T](top). The dependence 
of Ti on h is not monotonous: the shortest excitation time, 0.6 My, 
occurs at h = 0.75. Increasing asphericity, we reach » 2 My at 
h = 0.3. The other extreme would he h = I, but we stop at h = 0.99, 
because the results would be meaningless for a sphere. After cross- 
ing the separatrix, the angular momentum vector is driven towards 




Figure 1. Logarithmic plots of wobble damping/excitation times for ellip- 
soids with hi = /i2 = h. Top: time spent in the LAM with 9[ increasing from 
5° to 85°; bottom: time spent in the SAM with 83 decreasing from 85° to 
5° . Physical parameters - see the text. 



and we computed the damping times 7^3 taken to evolve from 
9!^ = 85° to 5°. This time (Fig.[T] bottom), damping times are much 
longer than in the LAM for the same shape. is as high as 258 My 
at A = 0.3 and systematically decreases to 6.5 My at h = 0.99. Thus, 
the total damping time consists mostly in and - fixing the semi- 
axis a - we find that triaxiality inhibits the total damping process, 
save for a quick passage through the LAM, where we find a mini- 
mum at h\ = 1-12 = 0.75. 

Integration limits in the above example have been wider than 
usually. For the reference with earlier works we have also computed 
sha pe factors Ds for the limits 45° < 0i < 85°, and 45° > 63 > 5°, 
like ISharma et al.1 ( l2005h . Fi gure |2] confronts Ds for a family of 
h\=h2 = h ellipsoids with spheroids having an appropriate ratio hi 
(SAM) or h2 (LAM) equal to 1, and the other one set as h. We note 
a systematic increase of Dj with increasing h for ellipsoids,similar 
to 73, but unlike Tj from Fig.[T] For oblate spheroids there is shal- 
low minimum of D3 close to hx: 0.9. These discrepancies are not 
essential and result from an arbitrary choice of integration limits 
that cut off some more or less prominent (depending on the shape) 
parts of an integrand. 

Integrated damping/excitation time is important for qualitative 
considerations, but if the joint action of inelastic dissipation and 
other torques is to be studied, the shape of ^sWs) becomes more 
interesting. Figures[3]and|4]demonstrate, how the ellipsoid's shapes 
affect ^1(6*1) and '^■^{6^). Each curve is normalized, i.e. "P., values 
are divided by w - the mean value of with respect to Os on the 
interval [0,;r/2]. The situation is a bit different in the LAM (Fig.O 
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0.3 0.4 0.5 0.6 0.7 

h 



0.8 0.9 



Figure 2. Logaiithmic plots of shape functions D,. Solid line: ellipsoids 
with h\ = h2 = ft; dashed line: spheroids with hi = h (LAM - top), or with 
h2=h (SAM - bottom). 




0, [deg] 

Figure 3. Normalized LAM energy dissipation rate functions 'Pi(6i)/w. 
Solid line: (fti.fti.w) = (0.7,1,3000); dashed: (0.7,0.7,4000); dotted: 
(0.3,0.7,9.6x 10'); dot-dashed: (0.7,0.3,5200). 



and SAM (Fig.Q. Departure from a spheroid (solid line) weakens 
the dissipation close to the principal axis and amplifies it for higher 
f^3 in the SAM. In the LAM, this is not the rule, as seen for hi = 0.3 
and h2 = 0.7. 




83 [deg] 

Figure 4. Normalized SAM energy dissipation rate functions ^,(9s)/w. 
Solid line: (fti,ft2,w) = (0.7,1,210); dashed: (0.7,0.7,530); dotted: 
(0.3,0.7, 1.6 X lO'*); dot-dashed: (0.7,0.3, 140). 



6 REDUCTION TO SPHEROID AND COMPARISON 
WITH OTHER WORKS 

In the previous section we have presented some results for 
spheroids. They could be computed by assuming hi or h2 suffi- 
ciently close to 1, but in the strict limit the expressions involve sin- 
gular factors. However this singularity is only apparent. The point 
is that although ki = qi =0 for h2 = 1, and % = 93 = for hi = 1, 
(hence all Pk{ks) = 0), but after we expand Pkiks) in powers of ks, 
substitute l |93l l and multiply by M,y, some factors ( 1 - /jj ) or ( 1 - /jj) 
cancel, leaving a well defined limit for a spheroid. The resulting 
SAM expression for hi = 1, /12 = h is 

8(l-/z2)sin2e3cose3 
T3(/j,03) = — ^ (2h^Ccos^9^+Ssm^ei),{iQ3) 



1,5 [\ + h^) 



where 



26 + 35/j2 

T3 + 20/J2' 

25 + 2Qlh^ + \6h^ 



(104) 



(105) 



: h, /i2 = 1> is simply 



l5 + Wh^ + 'ih'^ ' 
and the LAM expression for hi 

'¥i(h,6i) = -h'^'¥i{h-\0i). (106) 

The factor marks the difference between the LAM of the b = c < 
a spheroid (present work) and a prolate a = h < c spheroid used by 
other authors. 

We have found interesting to compare our solution with other 
published results. Remarkably, the latter can be reduced to same 
general form ( 11031 ) of 1'3, differing only with the particular expres- 
sions of the coefficients C and S . 

Let us begin with lEfroimskv & LazariaJ ( l200Cl() . Their solu- 
tion for a rectangular prism with semi-edges a = b and c = ah has 
the form of equation l ll03t with 



C 



1323 
128 ' 



105 



(107) 



Obviously, the prism has a higher volume than an ellipsoid with 
the same a and h. In particular, the volume integral jx^dV for a 
spheroid is smaller by a factor n/ 14 « 0.224. Thus, we propose to 
use 
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40 



e, [deg] 



Figure 5. "PsCSs) for a spheroid with h = c/a = 0.9. So l id line : present 
solution, dashed: v olume-scaled lEfroimsk v & Laz arianI j2bo(lt). dotted: 
iMolina et all j2003l) : dot-dashed: re-derived ISharma et al.l j2005l) for the 
equivalent Q definition. 




7T 1323 
C= — 



S, [deg] 

Figure 6. Same as Fig.|5]for h = 0.3. 
105 

^ " 14 16 ' 



(108) 



14 128 ' 

in equation il03t to make the comparison with a spheroid more 



iMolina et aD (l2003h obtained for a spheroid 



C= 1, 



2 



(109) 



Since our solutions differ only by the choice of boundary condi- 
tions, we adopt for comparison the values from equation ( 11091 ) with 

V = 1/4, i.e. C = 1 and S = 8/5^ 

The complete solution of ISharma et 'Zl ( I2005I) is known only 
indirectly, through the coefficients provided by Sharma (pri vate 
communication) and published in IVokrouhlickv et al] ( |2007[) . In 
these circumstances, we hav e derived C and 5 using the stress 
tensor from Appendix D of ISharma et al. I ( l2005b and following 
the recipe from their Sections 5 and 6 (i.e. the mainline solution, 
not the alternative model ). What we have obtained, agrees with 
IVokrouhlickv et alj ( l2007h . and reads 



C 
S 



13-)-20/!2 
25 + 20h^ + 1611"^ 



(110) 



where the factors in square brackets are the same as in our present 
solution JlOSt . The difference in denominators 2 and 4 is easily un- 
derstandable. First, the definition of Q used bv lSharma et ai]( l2005l) 
is based upon the mean value of energy, hence their energy dissipa- 
tion rate is twice as small as the one based upon our m ore common 
equation j4U . On the other hand. lSharma et all ([2005') do not apply 
the multiplier p in the sum given by our formula ( 147 1 . whereas S is 
directly related with the second harmonic of elastic energy. In these 
circumstances, we will use for the comparison 



C 



26-h35/i2 




13-i-20;i2 


■ -3 



25 + 20h^ + I6h'^ 



(111) 



i.e. the same C as in our solution dlOSt and a half of our S . 

Figures |5] and |6] present ^'^(O^) of the four solutions 
for oblate spheroids with h = c/a = 0.9 and 0.3. When the 
oblateness is moderate (Fig. O, there is a reasonable proxim- 
ity between the prese nt mo del and a volume- scal ed model of 
lEfroimskv & Lazaria nI j2000l). w hereas the models of lMolina et al.l 
I2OO3') and Sharma et al. I ( I2OO5I) have maxima lower by 25% than 
the present model and shifted with respect to each other by about 
15°. Curiously, increasing the oblateness (Fig. [6]l, we find a good 
agreement in the shape of the curves and dif ferences of the max- 
ima below 10%, with the notable exception o f ISharma et al] j2005t) 
that dissipates energy twice as slow as th e remaining models. A 
better agreement with IMolina et al.l ( |2003|) for smaller h is under- 
standable: the only difference with the present solution is due to 
different boundary conditions. They postulate a stress-free surface 
instead of the usual traction-free setup, and the deviatoric part of 
the stress tensor on the bounda ry in our solution decreases with 
h. Concerning the solution of S harma et al. U200l . we find a sys- 
tematic underestimation of 1^3 due to the missing multiplier of the 
second mode, with the ratio close to 2. 

Let us remind the claims of Shar ma et alj ( l200^ , that their so- 
lution gives damping times longer than other ones by factor 10 or 
more, and that the difference is due to incomplete or incorrect solu- 
tion of the elasticity problem in other papers. Why the differences 
in Figs. [5] and |6| are less drastic ? At this point, we feel obliged to 
observe that for any h, numerical v alues of shape factor s and damp- 
ing times published and plotted in ISharma et all ( l2005l Fig. 2) dif- 
fer from the ones resulting from his energy dissipation formulae 
by a constant factor n. We have found no trace of this discrepancy 
in the published equations, so the difference, most likely, should 
be attributed to a purely computational error. Together with the in- 
compatible definitmn_ofqu^ity factor Q, it means that all damping 
times from lSharma et al] ( |2005| ) should be divided by 2n, so they 
are no longer to be considered unusually high. On the other hand, 
damping times shorter than Burns-Safronov estimates claimed by 
lEfroimskv & LazarianI t2000) result partially (factor 1 4/;r) from us- 
ing an object (a prism) with a higher volume than any solid of rev- 
olution with the same ratio of axes. 



7 CONCLUSIONS 

Using the ensemble of standard assumptions, we have derived 
the stress tensor inside a freely rotating and self-gravitating 
ellipsoid. Writing abou t kno wn solutions to this problem, 
IWashabaugh & ScheeresI j2002h put meaningful quotation marks 
around the word available. Our present solution, given in the Ap- 
pendix|A]has a form which is probably compact enough to suppress 
the marks, especially for the principal axis rotation. Interestingly, 
the presented form of T does not involve singularity at the incom- 
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pressible limit v= where lWashabaugh & Scheerej ( l2002h had to 
use V = 0.499. 

The stress tensor has served us as t he basis for the en ergy 
dissipation model built along the lines that lEfroimskvl j2002h pro- 
posed, but left unaccomplished. However, the use of an ellipsoid 
instead of an Efroimsky-Lazarian rectangular prism permitted to 
avoid all objections related with partially satisfied boundary condi- 
tions and/or missing compatibility conditions. We have also found 
no reas ons to impose superfi cial conditions of a stress-free bound- 
ary like lMolina et alj j2003h . 

The solution hinges upon the use of Pki^l) series the Jacobi 
nome q. Their convergence is very good and even taking p < 4 in 
equations I l85l88t . guarantees at least three significant digits in the 
area under ^,(0^) for < 0j < n/2. Of course, the series are not le- 
gitimate exactly at q^ = 1 (i.e 9s = 90°). A n in-depth discussion of 
this limit was given bv lEfroimskvlj200lh . On the other hand, this 
state is not to be considered seriously, since any additional torque 
will trigger the emergence of a chaotic zone in the vicinity of sep- 
aratrices. 

In contrast to the results o f iMoUna et all ( |2003|) . we find that 
the role of compressibility in the energy dissipation process is 
marginal in the range of Poisson's ra tio ^ v < 4- . Appa rently, 
the stronger dependence obtained by iMolina et al.l ^20031) for a 
spheroid resulted from too strong boundary conditions. 

Investigating the spheroid as a particular case of our model, 
we have succeeded to resolve a major part of controversies con- 
cerning short damping times of lEfroimskv & LazarianI Ewdt) and 
very long ones according to ISharma et alj ( l2005h . In our opin- 
ion, the excess of energy dissipation rate over mainstream mod- 
els isjnostly_du£_to_aJimhe^^ of the body shape assumed 
by lEfroimskv & LazarianI ( l2000h . The shape factors reported by 
ISharma et al.l feOOSi) are overestimated mostly by an incompati- 
ble quality factor definition, a spurious factor n in their compu- 
tations, and a missing second mode multiplier. We find the objec- 
tions against bounda ry conditions used o r co mpatibility co n dition s 
violation, raised by iMolina et al.l ( l2003h or ISharma et"^ j2005l) . 
formally justified, yet we note that they affect the accuracy of the 
damping/excitation times by at most 50%. 
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APPENDIX A: ELLIPSOID STRESS TENSOR 

In order to shorten the expressions of the elements of A and A'-' in 
equation i32\ , we first introduce 



S'V* = _ ((- D'Bi 1 + (- iyS22 + (- D'^Bas) , 
and 

hi2 =hih2. 

Then, we can explicitly define 



^22 



(Al) 



(A2) 



(A3) 



© 2012 RAS, MNRAS 000.[Tlfl4l 



14 S. Breiter, A. Rozek and D. Vokrouhlicky 



= 3A23-/lf2523, 

= 2A„+/^72A22-/^7|A33-S""^ 

= -hl^Axx + /j|A22 + 2A33 - /Ji2S'°°, 

= 2Aii-;!i-2^22+V|A33-S°"', 

= 3Ai2-/2^Bi2, 

= -/i^Aii +2A22 + /l2^A33-/j^BlO°, 

= -hi{An+ hfA22 - h-jAii - fi'"") , 

= -h\^[lA2i-h\^B2i), 

= -hx[2An-h\^Bxi), 

= 2h\^[2h-'Ai2-hiBn), 

= 2h-,'{2h-,lA2i-hnB2i), 

= -hi[2h^^An-hi2Bi3), 

= -hn{2h-iUn-hiBi2}, 

= hihn{Au- hfA22 - h-^A33 + S'^O) , 

= -2/!72'^23 +^12^23, 

= -hi2{An- hfA22 + h-.lA^^ - sOio) , 

= 2h\[2h\\Aii-hnBn), 

= -hx2[2Ax2-h\Bx2). 



(A4) 
(A5) 
(A6) 
(A7) 
(A8) 
(A9) 
(AlO) 
(All) 
(A 12) 
(A 13) 
(A14) 
(A15) 
(A 16) 
(A17) 
(A 18) 
(A19) 
(A20) 
(A21) 
(A22) 
(A23) 



For the remaining 15 matrix elements that are not given above, we 
have (repeated index marks a pattern, no summation implied) 



A" - 



41 = 0. 



(A24) 

The off-diagonal 'central stress' elements are given directiy by 



A12 = 



1+y 




h\Bn 


2h'^^2hl + (3 + hl+hl2){l + v)^ 


2 ' 


(l + y)hl ] 


h 




2h] + (l + 3hl+hl2)(l+y) ) 




2 ' 






/?12S23 


2 + {h\ + h\^ + 3h\h\^){\ + v)_ 


2 



-423 



The rest is obtained from equations J35b with 



(A25) 



(A26) 



(A27) 



L„ = 


= -2-h\-h\-h\,^{\-h\)v, 


(A28) 


L12 = 


-- -\-h\-2h\ + h\2{\-h\)v. 


(A29) 


Ll3 = 


-- \+h\+h\ + h\{\ + 2h\ + 2h\2)v. 


(A30) 


L21 -- 


-- h\{l+h\ + h\) + {2 + 2h\+h\^)v. 


(A31) 


L22 = 


-- -h\(2 + h\ + h\)-{l-h\)v. 


(A32) 


^23 = 


-- -h\(l + hl + 2h\) + (\-h\)v. 


(A33) 


L31 = 


-- -2-h\2-h\2-h\{l-h\2)v, 


(A34) 


i32 = 


= 1 + /! 1 2 + + h\ (2 + h\ + 2/i j2 ) V, 


(A35) 


L33 = 


= -\-h\2-2h\2+h\{\-h\2)v. 


(A36) 


and 







\-h\v 



\-h\^y 



-hUl^h^^y ' 
h\^[h\-y) 



(A37) 



APPENDIX B: COEFFICIENTS M 

After setting the Poisson's ratio y=\, and defining 



N 



32 
35 



"12 



(Bl) 



where h\2 = h\h2, we obtain a compact form of three coefficients 
required in equations J82b and J83b 



Mi3 



M23 



M12 



w(l-/4)(l-/4) 
(l-/4)(l-/,4) 



2- 



5/zf 



5h\hl 



N 



8 + 5/!2 + 5/i22(l + 3/j2)^ 
5 



l5 + 5hl + hj^ (5 + 8^2)- 
The expression of the fourth one is more involved: 



(B2) 



(B3) 



(B4) 



Mo = 


3h\N,jri ' ' 


(B5) 


Using ; 


m auxiliary variable 








(B6) 


we can 


compress Nj to read 




No = 


225(^-1), 


(B7) 


Ni = 


6(l + /;^)(29^-21), 


(B8) 


N2 = 


(31^^ + 82^ -62), 


(B9) 


N2 = 


h\[l+hfj [-92 f + 305^ - 2 16) , 


(BIO) 


N4 = 


/("f (3 1^^-34 1^2 + 99^ + 295), 


(BU) 


Ns = 


h\{l + /ii)(l74f^ - 1012^2 + 1 185^-458) , 


(B12) 


A'6 = 


h\ (225^ - + 2412^2 - 1409^ - 124) , 


(B13) 


yVv = 


h\{l + hfj[225f - 1 179^2 + 1376^-368) , 


(B14) 


A'8 = 


h\ (3^ - 4) {l5f- - 292^ + 64) , 


(B15) 


Ni) = 


4S^-57 + h]hl[4Sf-l\9^+ lOO) 
+hl[\ + h\) (32^ - 23 + h^hl (39^ - 44)) 






+ 16/i|;z^(3^-4). 


(B16) 
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